##############################################################################
##############################################################################

###              R code to replicate findings presented in                 ###
###                 "The Index of Emancipative Values:                     ### 
###                  Measurement Model Misspefication"                     ###    
###        (forthcoming in the American Political Science Review)          ### 
                          
##############################################################################
##############################################################################



##############################################################################

###                   I. DATA AND REQUIRED R PACKAGES                      ###

##############################################################################

library(foreign)
library(semPlot)
library(stargazer)
library(MplusAutomation)
library(car)
library(lmtest)


setwd( "C:/Documents and Settings/User")
#WVS<- read.spss("WVS_Longitudinal_1981_2014_spss_v2015_04_18.sav",
#                 to.data.frame=T, 
#                 use.value.labels=T)

WVS_3_6 <- read.spss("wvs_evi_3_6.sav",
                  to.data.frame=T, 
                  use.value.labels=T)

#WVS_3_6 <- subset(WVS, 
#WVS$S002 == 3|WVS$S002 == 4|WVS$S002 == 5|WVS$S002 == 6, 
#select=c(S002, S003, S017, S019,
#C001, D059, D060,
#A029, A034, A042,
#F118, F120, F121,
#E001, EOO2, E003, EOO4))

## Relevant variable labels
## S002: The wave
## S003: The country
## S017: The original weight
## S019: The original weight equilibrated to an homogeneous N of 1500 for each 
## 
## Equality: C001, D059, D060  (job, leadership, education)
## 
## Autonomy: A029, A034, A042 (Autonomy, Imagination, Obedience)
## 
## Freedom: F118, F120, F121  (Homosexuality, Abortion, Divorce)
## 
## Voice: 
## E001 3 (1st) and EOO2 3 (2nd): Aims of Country: more say (local)
## E003 2 (1st) and EOO4 2 (2nd): Aims of Respondent: more say (national)
## E003 4 (1st) and EOO4 4 (2nd): Aims of Respondent: freedom of speech 

   


##############################################################################

###                  II. SURVEY ITEMS MEASURING THE EVI                    ###

##############################################################################

### AUTONOMY ###

WVS_3_6$Independence <- recode(as.numeric(WVS_3_6$A029), "1 = 0; 2=1")
WVS_3_6$Imagination <- recode(as.numeric(WVS_3_6$A034), "1 = 0; 2=1")
WVS_3_6$Obedience <- recode(as.numeric(WVS_3_6$A042), "1 = 1; 2=0")

### EQUALITY ###

WVS_3_6$Jobs <-recode(as.numeric(WVS_3_6$C001), "1=0; 2=0.5; 3=1")
WVS_3_6$Leaders <-recode(as.numeric(WVS_3_6$D059), "1=0; 2=0.33; 3=0.75; 4=1")
WVS_3_6$Education <-recode(as.numeric(WVS_3_6$D060), "1=0; 2=0.25; 3=0.75; 4=1")

### FREEDOM ###

WVS_3_6$Homosexuality <-(as.numeric(WVS_3_6$F118)-1)/9
WVS_3_6$Abortion <-(as.numeric(WVS_3_6$F120)-1)/9
WVS_3_6$Divorce <-(as.numeric(WVS_3_6$F121)-1)/9

### VOICE ###
WVS_3_6$Say_nat[as.numeric(WVS_3_6$E003) != 2 & as.numeric(WVS_3_6$E004) != 2] <-0
WVS_3_6$Say_nat[as.numeric(WVS_3_6$E004) == 2] <-0.5
WVS_3_6$Say_nat[as.numeric(WVS_3_6$E003) == 2] <-1

WVS_3_6$Speech[as.numeric(WVS_3_6$E003) != 4 & as.numeric(WVS_3_6$E004) != 4] <-0
WVS_3_6$Speech[as.numeric(WVS_3_6$E004) == 4] <-0.5
WVS_3_6$Speech[as.numeric(WVS_3_6$E003) == 4] <-1

WVS_3_6$Say_local[as.numeric(WVS_3_6$E001) != 3 & as.numeric(WVS_3_6$E002) != 3] <- 0
WVS_3_6$Say_local[as.numeric(WVS_3_6$E002) == 3] <-0.5
WVS_3_6$Say_local[as.numeric(WVS_3_6$E001) == 3] <-1




##############################################################################

###             III. FIGURE 1: MEASUREMENT MODEL FOR THE EVI               ###

##############################################################################


# define a CFA measurement model for the EVI
EVI <- 'Autonomy  =~ Independence + Imagination + Obedience 
              Equality =~ Jobs + Leaders + Education
              Choice  =~ Homosexuality + Abortion + Divorce 
              Voice =~ Speech + Say_nat + Say_local 
              EVI =~ Autonomy + Equality + Choice + Voice'

# plot the measurement model

tiff('Figure 1.tiff', units = "cm", width = 15.5, height = 9, res = 1000)

semPaths(semPlotModel(EVI), curvePivot = TRUE,
residuals = FALSE,
nodeLabels= c("Inde-\npendence", "Imagi-\nnation", "Obedience", # set node labels
"Jobs", "Leaders", "Education",
"Homo-\nsexuality", "Abortion", "Divorce",
"Freedom \nof speech", "More Say:\n National", "More Say:\n Local", 
"Autonomy", "Equality", "Choice", "Voice", "Index of\nEmancipative\nValues" ),
optimizeLatRes = TRUE,
sizeMan = 9,
sizeMan2 = 5,
sizeLat = 12,
mar = c(5,1,7,1),
label.cex= 0.6, 
label.scale= FALSE, # label size is not matched to the size of node
edge.width= 1.5,  
edge.color= "black")  # color of edge lines


dev.off()



##############################################################################

###               IV.  DATA PREPARATION FOR MPLUS ANALYSIS                 ###
###       (Table 1; Tables A1 and B1-B4 in Supplementary Materials)        ###   

##############################################################################


### Create a new variable "Cultural zone" (as defined by Welzel [2013]) 

WVS_3_6$Zone[WVS_3_6$S003 == "Algeria" | 
WVS_3_6$S003 == "Iran" | 
WVS_3_6$S003 == "Iraq" |
WVS_3_6$S003 == "Libya" | 
WVS_3_6$S003 == "Saudi Arabia" |
WVS_3_6$S003 ==  "Egypt" |
WVS_3_6$S003 ==  "Jordan" | 
WVS_3_6$S003 == "Morocco" |
WVS_3_6$S003 == "Turkey"|
WVS_3_6$S003 == "Palestine" |
WVS_3_6$S003 == "Libya"|
WVS_3_6$S003 ==  "Yemen" |
WVS_3_6$S003 ==  "Lebanon" |
WVS_3_6$S003 == "Tunisia" |
WVS_3_6$S003 == "Bahrain" |  
WVS_3_6$S003 == "Qatar"|            
WVS_3_6$S003 == "Kuwait"] <- "Islamic East"

WVS_3_6$Zone[WVS_3_6$S003 == "Bangladesh" | 
WVS_3_6$S003 =="Indonesia"|
WVS_3_6$S003 =="Pakistan" | 
WVS_3_6$S003 =="India"|
WVS_3_6$S003 =="Malaysia" |
WVS_3_6$S003 == "Philippines"|
WVS_3_6$S003 =="Singapore"|
WVS_3_6$S003 =="Thailand"] <- "Indic East"

WVS_3_6$Zone[WVS_3_6$S003 == "Viet Nam"|
WVS_3_6$S003 =="China"|
WVS_3_6$S003 =="Hong Kong"|
WVS_3_6$S003 =="South Korea"|
WVS_3_6$S003 =="Taiwan"|
WVS_3_6$S003 =="Japan"] <- "Sinic East"  

WVS_3_6$Zone[WVS_3_6$S003 == "Azerbaijan"|
WVS_3_6$S003 =="Albania"|
WVS_3_6$S003 =="Armenia"|
WVS_3_6$S003 =="Belarus"|
WVS_3_6$S003 =="Bosnia"|
WVS_3_6$S003 =="Georgia"|
WVS_3_6$S003 =="Kyrgyzstan"|
WVS_3_6$S003 =="Russia"|
WVS_3_6$S003 =="Serbia"|
WVS_3_6$S003 =="Kazakhstan"|
WVS_3_6$S003 =="Uzbekistan"|
WVS_3_6$S003 =="Ukraine"|
WVS_3_6$S003 =="Macedonia"|
WVS_3_6$S003 =="Moldova"|
WVS_3_6$S003 =="Romania"|
WVS_3_6$S003 =="Bulgaria"|
WVS_3_6$S003 =="Serbia and Montenegro"|
WVS_3_6$S003 =="Montenegro"] <- "Orthodox East" 

WVS_3_6$Zone[WVS_3_6$S003 == "Cyprus"|
WVS_3_6$S003 =="Greece"|
WVS_3_6$S003 =="Israel"|
WVS_3_6$S003 =="Andorra"|
WVS_3_6$S003 =="Austria"|
WVS_3_6$S003 =="Belgium"|
WVS_3_6$S003 =="France"|
WVS_3_6$S003 =="Ireland"|
WVS_3_6$S003 =="Italy"|
WVS_3_6$S003 =="Luxemburg"|
WVS_3_6$S003 =="Malta"|
WVS_3_6$S003 =="Portugal"|
WVS_3_6$S003 =="Spain"] <- "Old West" 

WVS_3_6$Zone[WVS_3_6$S003 == "Denmark"|
WVS_3_6$S003 =="Finland"|
WVS_3_6$S003 =="Germany"|
WVS_3_6$S003 =="Iceland"|
WVS_3_6$S003 =="Netherlands"|
WVS_3_6$S003 =="Norway"|
WVS_3_6$S003 =="Sweden"|
WVS_3_6$S003 =="Switzerland"|
WVS_3_6$S003 =="Great Britain"] <- "Reformed West" 

WVS_3_6$Zone[WVS_3_6$S003 == "Australia"|
WVS_3_6$S003 =="Canada"|
WVS_3_6$S003 =="New Zealand"|
WVS_3_6$S003 =="United States"] <- "New West" 

WVS_3_6$Zone[WVS_3_6$S003 == "Croatia" |
WVS_3_6$S003 =="Latvia"|
WVS_3_6$S003 =="Lithuania"|
WVS_3_6$S003 =="Estonia"|
WVS_3_6$S003 =="Czech Rep."|
WVS_3_6$S003 =="Hungary"|
WVS_3_6$S003 =="Poland"|
WVS_3_6$S003 =="Slovakia"|
WVS_3_6$S003 =="Slovenia"] <- "Returned West"  

WVS_3_6$Zone[WVS_3_6$S003 == "Guatemala"|
WVS_3_6$S003 =="Venezuela"|
WVS_3_6$S003 =="Ecuador"|
WVS_3_6$S003 =="Colombia"|
WVS_3_6$S003 =="Mexico"|
WVS_3_6$S003 =="Peru"|
WVS_3_6$S003 =="Brazil"|
WVS_3_6$S003 =="Chile"|
WVS_3_6$S003 == "Puerto Rico"|
WVS_3_6$S003 =="Dominican Rep."|
WVS_3_6$S003 =="El Salvador"|
WVS_3_6$S003 =="Argentina"|
WVS_3_6$S003 =="Trinidad and Tobago"|
WVS_3_6$S003 =="Uruguay"] <- "Latin America" 

WVS_3_6$Zone[WVS_3_6$S003 == "Burkina Faso"|
WVS_3_6$S003 =="Ghana"|
WVS_3_6$S003 =="Ethiopia"|
WVS_3_6$S003 =="Nigeria"|
WVS_3_6$S003 =="Rwanda"|
WVS_3_6$S003 =="Tanzania"|
WVS_3_6$S003 =="Uganda"|
WVS_3_6$S003 =="Zimbabwe"|
WVS_3_6$S003 =="Mali"|
WVS_3_6$S003 =="Zambia"|
WVS_3_6$S003 =="South Africa"] <- "sub-Saharan Africa"  


WVS_3_6$Zone <- as.factor(WVS_3_6$Zone)
WVS_3_6$S019 <- as.numeric(WVS_3_6$S019) # Country weights


# Define which indicators are ordered categorical variable
# Note that MPLUS doesn't treat a variable as categorical
# if it contains non-integer values
WVS_3_6$Independence_Cat <- as.factor(WVS_3_6$Independence)
WVS_3_6$Imagination_Cat <- as.factor(WVS_3_6$Imagination)
WVS_3_6$Obedience_Cat <- as.factor(WVS_3_6$Obedience)
WVS_3_6$Jobs_Cat <- as.ordered(WVS_3_6$Jobs)
WVS_3_6$Leaders_Cat <- as.ordered(WVS_3_6$Leaders)
WVS_3_6$Education_Cat <- as.ordered(WVS_3_6$Education)
WVS_3_6$Speech_Cat <- as.ordered(WVS_3_6$Speech)
WVS_3_6$Say_nat_Cat <- as.ordered(WVS_3_6$Say_nat)
WVS_3_6$Say_local_Cat <- as.ordered(WVS_3_6$Say_local)


# Create an MPLUS input data file

prepareMplusData(WVS_3_6, "wvs3_6.dat",
keepCols = c("S002", "S003", "S019", "Zone", 
"Independence_Cat", "Imagination_Cat", "Obedience_Cat", 
"Jobs_Cat", "Leaders_Cat", "Education_Cat",
"Homosexuality","Abortion", "Divorce", 
"Speech_Cat", "Say_nat_Cat", "Say_local_Cat"))


############################################################################

###                        III. TABLES 2 AND 3                           ###  
###           (Associations between the EVI, Its Components,             ### 
###               and the Index of Effective Democracy)                  ### 

############################################################################


### Select data for the WVS 6th round 

WVS_6 <- WVS_3_6[WVS_3_6$S002 == "2010-2014",]


### Compute individual scores on the EVI and each of its particular component

# Autonomy
WVS_6$aut <- rowMeans(data.frame(WVS_6$Independence, 
WVS_6$Imagination, WVS_6$Obedience), na.rm = T)

# Equality
WVS_6$eq <- rowMeans(data.frame(WVS_6$Jobs,WVS_6$Leaders, 
WVS_6$Education), na.rm = T)

# Choice
WVS_6$ch<- rowMeans(data.frame(WVS_6$Homosexuality, 
WVS_6$Abortion, WVS_6$Divorce), na.rm = T)

# Voice
WVS_6$vc <- rowMeans(data.frame(WVS_6$Speech, 
WVS_6$Say_nat, WVS_6$Say_local), na.rm = T)

# The whole index
WVS_6$emval<-rowMeans(data.frame(WVS_6$aut, WVS_6$eq, WVS_6$ch, WVS_6$vc), na.rm =T)



###  Compute country mean scores on the EVI and each of its particular component

aut.means <- aggregate(WVS_6$aut, by = list(WVS_6$S003), FUN = mean, na.rm = T)
eq.means <- aggregate(WVS_6$eq, by = list(WVS_6$S003), FUN = mean, na.rm = T)
ch.means <- aggregate(WVS_6$ch, by = list(WVS_6$S003), FUN = mean, na.rm = T)
vc.means <- aggregate(WVS_6$vc, by = list(WVS_6$S003), FUN = mean, na.rm = T)
EVI.means <- aggregate(WVS_6$emval, by = list(WVS_6$S003), FUN = mean, na.rm = T)


# Estimated [standaridized] latent mean scores on "choice" 
# for 58 countries covered by the 6th wave of the WVS 
# (approximate metric and scalar invariance is assumed for each inidcator)
# (see folder "MPLUS Pro-Choice Input", MPLUS input file "Wave 6")

latent.means <- c(-0.290, -0.743, 1.022, 
1.218, -1.238, -1.090, 
0.742, 0.391, 1.128, 
-0.483, 0.197, -0.199, 
0.801, -0.454, 0.717,
-1.016, -1.066, 0.795, 
-1.771, -0.051, -0.213,
-0.529, 0.671, 0.063, 
-1.606, 0.067,  NA,
-0.951, -0.402, -1.090,
-0.977, 0.085, -0.390, 
1.145, 0.967, -1.307, 
-2.742, 0.061, -0.210, 
0.198, -1.062, -0.275, 
0.718, -2.104, -0.121, 
1.119, 0.068, -0.982, 
1.330, 1.914, -0.670, 
-0.419, -0.635, -0.669, 
0.380, NA, 0.806, 
1.150, -0.918, -0.542)

# *** NAs are Kuwait and Egypt. In these countries, 
#  one or more items measuring pro-choice orientations 
#  were not included in the national WVS questionnaires and 
#  therefore latent means for �choice� could not be estimated)



### Data to compute the Index of Effective Democracy ###

# Country Labels
countries <- c("Algeria", "Azerbaijan", "Argentina", 
"Australia", "Bahrain", "Armenia",
"Brazil", "Belarus", "Chile", 
"China", "Taiwan", "Colombia", 
"Cyprus", "Ecuador", "Estonia", 
"Georgia", "Palestine", "Germany",
"Ghana", "Hong Kong", "India",
"Iraq", "Japan", "Kazakhstan",
"Jordan", "South Korea", "Kuwait",
"Kyrgyzstan", "Lebanon", "Libya", 
"Malaysia", "Mexico", "Morocco", 
"Netherlands", "New Zealand", "Nigeria",
"Pakistan", "Peru", "Philippines", 
"Poland", "Qatar", "Romania",
"Russia", "Rwanda", "Singapore", 
"Slovenia", "South Africa", "Zimbabwe",           
"Spain", "Sweden", "Thailand",
"Trinidad and Tobago", "Tunisia", "Turkey", 
"Ukraine", "Egypt", "United States", 
"Uruguay", "Uzbekistan", "Yemen")


# 2014 Rule of Law Score for 60 Countries #
# Source: http://databank.worldbank.org/data/reports.aspx?source=worldwide-governance-indicators#

rule.of.law <- c(-0.77, -0.67, -0.89, 
1.92, 0.41, -0.37,
-0.05, -0.84, 1.43,
-0.41, 1.19, -0.29,
1.08, -1.06, 1.37,
0.19, -0.39, 1.86,
0.05, 1.86, -0.06,
-1.33, 1.60, -0.60, 
0.46, 0.99, 0.02, 
-0.93, -0.77, -1.51, 
0.59, -0.42, -0.07, 
1.98, 2.00, -1.05, 
-0.76, -0.52, -0.32, 
0.84, 0.86, 0.17, 
-0.74, 0.06, 1.82, 
1.00, 0.18, -1.43, 
0.95, 1.99, -0.19, 
-0.09,-0.11, 0.00, 
-0.79, -0.66, 1.61,
0.71, -1.13, -1.20)


# 2014 Control of Corruption Score for 60 Countries #
# Source: http://databank.worldbank.org/data/reports.aspx?source=worldwide-governance-indicators#

control.of.corruption <- c(-0.60, -1.02, -0.54,
1.85, 0.28, -0.52, 
-0.34, -0.30, 1.49,
-0.34, 0.81, -0.37,
1.08, -0.75, 1.30,
0.79, -0.28, 1.84,
-0.19, 1.62, -0.43,
-1.33, 1.69, -0.83,
0.14, 0.55, -0.24, 
-1.13, -1.04, -1.56,
0.41, -0.76, -0.27, 
1.99, 2.25, -1.27, 
-0.83, -0.56, -0.44, 
0.64, 0.99, -0.11,
-0.92, 0.76, 2.07, 
0.73, -0.06, -1.39, 
0.63, 2.15, -0.45, 
-0.41, -0.04, -0.15, 
-0.99, -0.62, 1.38, 
1.38, -1.19, -1.56)


# 2014 Freedom House Scores for 60 Countries #
# Source: https://freedomhouse.org/sites/default/files/FIW2014%20Booklet.pdf
freedoms <- matrix(c(6, 5, 6, 6, 2, 2,
1, 1, 6, 6, 5, 4,
2, 2, 7, 6, 1, 1,
7, 6, 1, 2, 3, 4,
1, 1, 3, 3, 1, 1,
3, 3, 6, 5, 1, 1,
1, 2, 5, 2, 2, 3,
5, 6, 1, 1, 6, 5,
6, 5, 2, 2, 5, 5, 
5, 5, 5, 4, 4, 5, 
4, 4, 3, 3, 5, 4, 
1, 1, 1, 1, 4, 4, 
4, 5, 2, 3, 3, 3, 
1, 1, 6, 5, 2, 2, 
6, 5, 6, 5, 4, 4, 
1, 1, 2, 2, 5, 6, 
1, 1, 1, 1, 4, 4, 
2, 2, 3, 3, 3, 4, 
4, 3, 6, 5, 1, 1, 
1, 1, 7, 7, 6, 6), ncol = 2, byrow= TRUE)


### Construct Democratic Rights Index (Alexander et al. 2012) ###
DRI <-  (14 - (freedoms[,1] + freedoms[,2]))/0.12
COS <- (rule.of.law + control.of.corruption)/200
RLI <- (COS - min(COS))/(max(COS) - min(COS))
EDI <- (DRI*RLI)/100


###       Table II: Pairwise Correlations between Effective       ###
###    Democracy and Various Components of Emancipative Values    ###

round(cor(data.frame(EDI, DRI, latent.means, ch.means[,2],  
aut.means[,2], vc.means[,2], eq.means[,2], EVI.means[,2]),
use = "complete.obs")[1,3:8],3)

round(cor(data.frame(EDI, DRI, latent.means, ch.means[,2],  
aut.means[,2], vc.means[,2], eq.means[,2], EVI.means[,2]), method = "spearman",
use = "complete.obs")[1,3:8],3)

round(cor(data.frame(EDI, DRI, latent.means, ch.means[,2],  
aut.means[,2], vc.means[,2], eq.means[,2], EVI.means[,2]),method = "kendall",
use = "complete.obs")[1,3:8],3)



###    Table III: Regressing Effective Democracy on    ###
###         Components of Emancipative Values          ###

model1 <- lm(EDI ~ vc.means[,2])
summary(model1)
ncvTest(model1)
bptest(model1)
outlierTest(model1)
coeftest(model1, vcov = hccm) 

model2 <- lm(EDI ~ aut.means[,2])
summary(model2)
ncvTest(model2)
bptest(model2)
outlierTest(model2)
coeftest(model2, vcov = hccm) 


model3 <- lm(EDI ~ vc.means[,2] + ch.means[,2])
summary(model3)
vif(model3)
outlierTest(model3)
ncvTest(model3)
bptest(model3)
coeftest(model3, vcov = hccm) 

model4 <- lm(EDI ~ aut.means[,2] + ch.means[,2])
summary(model4)
vif(model4)
outlierTest(model4)
ncvTest(model4)
bptest(model4)
coeftest(model4, vcov = hccm) 

model5 <- lm(EDI ~ vc.means[,2] + eq.means[,2])
summary(model5)
vif(model5)
outlierTest(model5)
ncvTest(model5)
bptest(model5)
coeftest(model5, vcov = hccm) 

model6 <- lm(EDI ~ aut.means[,2] + eq.means[,2])
summary(model6)
vif(model6)
outlierTest(model6)
ncvTest(model6)
bptest(model6)
coeftest(model6, vcov = hccm) 

model7 <- lm(EDI ~ ch.means[,2] + eq.means[,2])
summary(model7)
vif(model7)
outlierTest(model7)
ncvTest(model7)
bptest(model7)
coeftest(model7, vcov = hccm) 

model8 <- lm(EDI ~ EVI.means[,2])
summary(model8)
outlierTest(model8)
ncvTest(model8)
bptest(model8)
coeftest(model8, vcov = hccm) 


stargazer(model1, model2, model3, model4, model5, model6, model7, model8,
type = "html", title = "Table 3. Regressing Effective Democracy on Components of Emancipative Values",
out= "evi_validity.htm", style = "apsr",
covariate.labels = c("Voice", "Autonomy", "Choice", "Equality", "EVI"),
no.space = T, single.row = T,
star.cutoffs = c(0.05, 0.01, 0.005),
notes.align = "c")

